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Abstract 

Initial data are the starting point for any numerical simulation. In 
the case of numerical relativity, Einstein's equations constrain our choices 
of these initial data. We will examine several of the formalisms used for 
specifying Cauchy initial data in the 3 + 1 decomposition of Einstein's 
equations. We will then explore how these formalisms have been used in 
constructing initial data for spacetimes containing black holes and neutron 
stars. In the topics discussed, emphasis is placed on those issues that are 
important for obtaining astrophysically realistic initial data for compact 
binary coalescence. 
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1 Introduction 



The goal of numerical relativity is to study spacctimes that cannot be studied 
by analytic means. The focus is therefore primarily on dynamical systems. Nu- 
merical relativity has been applied in many areas: cosmological models, critical 
phenomena, perturbed black holes and neutron stars, and the coalescence of 
black holes and neutron stars, for example. In any of these cases, Einstein's 
equations can be formulated in several ways that allow us to evolve the dynam- 
ics. While Cauchy methods have received a majority of the attention, charac- 
teristic and Reggi calculus based methods have also been used. All of these 
methods begin with a snapshot of the gravitational fields on some hypersurface, 
the initial data, and evolve these data to neighboring hypersurfaces. 

The focus of this review is on the initial data needed for Cauchy evolutions 
of Einstein's equations. These initial data cannot be freely specified in their 
entirety. Rather they are subject to certain constraints which must be satisfied. 
Because of the nonlinearity of Einstein's equations, there is no unique way of 
choosing which pieces of the initial data can be freely specified and which are 
constrained. In §2, I will look at the various formalisms that exist for expressing 
the initial- value equations. I will use the 3 -fl (or ADM) decomposition of 
Einstein's equations throughout and I begin the review with a brief introduction 
of this in §2.1. In §2.2, §2.3, and §2.4, I will explore the most important and 
widely used decompositions that have been developed. 

In the remainder of the review, I will focus on initial data for black holes and 
neutron stars. Section 3 deals with black- hole initial data, which has received 
considerable attention over the years. The evolution of black-hole spacetimes is 
of particular importance because it allows for the study of pure geometrodynam- 
ics. The spacetime is either pure vacuum or any matter can be hidden behind 
an event horizon. Thus, black- hole spacetimes allow us to study the dynamics 
of gravity alone. In §3.1, I begin with a review of some of the classic analytic 
solutions of the initial-data equations. Section 3.2 covers current methods for 
constructing general multi-black-hole initial data sets and also explores some 
of the limitations of these data sets. Section 3.3 deals with a class of black- 
hole initial data that has received little attention until recently and which I feel 
may play an important role in future work on constructing black-hole initial 
data. Finally, §3.4 examines the issue of constructing black-hole initial data for 
quasiequilibrium binaries. 

This last topic is of extreme importance. Several laser interferometer gravi- 
tational wave detectors will become operational in the near future and the coa- 
lescence of black- hole binaries is considered to be one of the strongest candidates 
for detection by the earliest generation of detectors. The chances of detecting 
these events and then unraveling the information contained in the gravitational 
wave signals will be greatly increased if we have accurate numerical simulations 
of black-hole binary coalescence. While post-Newtonian techniques can be used 
to simulate the inspiral of a compact binary system, the final plunge and co- 
alescence must be simulated numerically. Astrophysically realistic initial data 
will be needed before these simulations can provide reliable results. 
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The final plunge and coalescence of neutron-star binaries must also be sim- 
ulated numerically. Section 4 examines the construction of neutron-star initial 
data. The major difference between black-hole and neutron-star initial data is 
the need to deal with the neutron star's matter. In this section, I will deal 
with the neutron star matter in general, without considering any particular 
equation of state. In §4.1, I look at the issue of hydrostatic equilibrium of mat- 
ter. Section 4.2 takes a brief look at constructing equilibrium initial data for 
isolated neutron stars. Finally, §4.3 examines the issues involved in construct- 
ing quasiequilibrium neutron-star binary initial data. In particular, we look at 
recent advances in constructing irrotational fluid models. 

This article is far from a complete review of all the work that has been done 
on initial data for numerical relativity. I offer my apologies to all whose work I 
have not included. In particular, I have not covered any of the issues associated 
with existence and uniqueness of solutions to the constraint equations, nor dealt 
extensively with the issue of asymptotic falloff rates. I welcome correspondence 
on topics that you feel should be covered in this review, references that you 
think should be included, and the inevitable typographical errors. 

1.1 Conventions 

I will use a 4- metric g^u with signature ( — h ++). Following the MTW [79] 
conventions, I define the Riemann tensor as 

The Ricci tensor is defined as 

R^v = R ixau-, (2) 

and the Einstein tensor as 

Gjxv = Rjxv 29ni'R' {^) 

A spatial 3-metric will be written as jij and the Riemann and Ricci tensors 
associated with it will be defined by (1) and (2). Four-dimensional indices will 
be denoted by Greek letters, while 3-dimensional indices will be denoted by 
Latin letters. 

To avoid confusion, the covariant derivative and Ricci tensor associated with 
jij will be written with over-bars — Vj and Rij . I will also frequently deal with 
an auxiliary 3-dimensional space with a metric that is conformally related to 
the metric 7ij of the physical space. The metric for this space will be denoted 
7jj. The covariant derivative and Ricci tensor associated with this metric will 
be written with tildes — Vj and Rij. Other quantities that have a conformal 
relationship to quantities in the physical space will also be written with a tilde 
over them. 
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2 The Initial- Value Equations 



Einstein's equations, G^v = S-kGT^^, represent ten independent equations. 
Since there are ten equations and ten independent components of the four- 
metric it seems that we have tlie same number of equations as imknowns. 
From the definition of the Einstein tensor (3), we see that these ten equations 
are linear in the second derivatives and quadratic in the first derivatives of the 
metric. We might expect that these ten second-order equations represent evolu- 
tion equations for the ten components of the metric. However, a close inspection 
of the equations reveals that only six of the ten involve second time-derivatives 
of the metric. The remaining four equations are not evolution equations. In- 
stead, they are constraint equations. The full system of equations is still well 
posed, however, because of the Bianchi identities 

V^G""' = 0. (4) 

The four constraint equations appear as a result of the general covariance of 
Einstein's theory, which gives us the freedom to apply general coordinate trans- 
formations to each of the four coordinates and leave the interval, 

ds2 = c,^^da;^da;'', (5) 

unchanged. 

If we consider Einstein's equations as a Cauchy problem^, we find that the 
ten equations separate into a set of four constraint or initial-value equations, 
and six evolution or dynamical equations. If the four initial- value equations are 
satisfied on some spacelike hypersurface, which we can label with t = 0, then 
the Bianchi identities (4) guarantee that the evolution equations preserve the 
constraints on neighboring spacelike hypersurfaces. 

2.1 Initial Data 

In the Cauchy formulation of Einstein's equations, we begin by foliating the 

4-dimensional manifold as a set of spacelike, 3-dimensional hypersurfaces (or 
slices) {S}. These slices are labeled by a parameter t or, more simply, each 
slice of the 4-dimensional manifold is a f = constant hypersurface. Following 
the standard 3 + 1 decomposition [6, 110, 37], we let n'^ be the future-pointing 
timelike unit normal to the slice, with 

n'^ = -aV^i. (6) 

Here, (y IS c'f lUed the lapse function (frequently denoted N in the literature). 
The scalar lapse function sets the proper interval measured by observers as they 
move between slices on a path that is normal to the hypersurface (so-called 
normal observers): 

^slalong = "dt (7) 

^We could also formulate Einstein's equations as a chara<;teristic initial- value problem, but 
we will not pursue that approach in this paper. 
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Of course, there is no reason that observers must move along a path normal to 
the hypersurface. In general, we can define the time vector as 

= an" + (3^, (8) 

where 

f3^n^ = 0. (9) 

Here, (3^ is called the shift vector (frequently denoted in the literature). 

Because of (9), has only three independent components and is a spatial 
vector, tangent to the hypersurface on which it resides. At this point, it is 
convenient to introduce a coordinate system adapted to the foliation {E}. Let 

be the spatial coordinates in the slice. The fourth coordinate, t, is the 
parameter labeling each slice. With this adapted coordinate system, wc find 
that 3-dimensional coordinate values remain constant as we move between slices 
along the direction (8). The four parameters, a and are a manifestation of 
the 4-dimensional coordinate invariance, or gauge freedom, in Einstein's theory. 
If we let 7y represent the metric of the spacelike hypersurfaces, then we can 
rewrite the interval (5) as 

ds^ = -a^dt^ + 7y (dx' + I3'dt){dx^ + jS^dt). (10) 

In the Cauchy formulation of Einstein's equations, jij is regarded as the 
fundamental variable and values for its components must be given as part of a 
well-posed initial-value problem. Since Einstein's equations are second order, 
we must also specify something like a time derivative of the metric. For this, 
we use the second fundamental form, or extrinsic curvature, of the slice, Kij , 
defined^ by 

K,, = -i/:„7,;,. (11) 

where £„ denotes the Lie derivative along the n'' direction. 

Together, jij and Kij are the minimal set of initial data that must be spec- 
ified for a Cauchy evolution of Einstein's equations. The metric 7ij on a hy- 
persurface is induced on that surface by the 4- metric g^i,. This means that the 
values jij receives depend on how E is embedded in the full spacetime. In order 
for the foliation of slices {S} to fit into the higher-dimensional space, they must 
satisfy the Gauss-Codazzi-Ricci conditions. Combining these conditions with 
Einstein's equations, and using (10), the six evolution equations become 

dtKij = a[R^J-2KuK'J+KK^J-8^TGS^J+i^TGJij{S-p)] (12) 

Here, Vj is the spatial covariant derivative compatible with 7ij , Rij is the Ricci 
tensor associated with 7^ , K = Kl, p is the matter energy density, 5^ is 

diflferent sign choice for defining the extrinsic curvature is sometimes found in the 
literature. 
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the matter stress tensor, and S = Sl-^ We have also used the fact that in 
our adapted coordinate system, = dt- The set of second-order evolution 
equations is completed by rewriting the definition of the extrinsic curvature 

(11) as 

daij = -2aK,j + + VjA. (13) 

Equations (12) and (13) are a first-order representation of a complete set of 
evolution equations for given initial data 7jj and Kij. However, the data cannot 

be freely specified in their entirety. The four constraint equations, following the 
same procedure outlined above for the evolution equations, become 

R + K^- KijK'^ = WnGp (14) 

and 

{K'^ - -f'^K) = SnGf. (15) 

Here, R= R\ and f is the matter momentum density. Equation (14) is referred 
to as the Hamiltonian or scalar constraint, while (15) are referred to as the 
momentum or vector constraints. Valid initial data for the evolution equations 

(12) and (13) must satisfy this set of constraints. And, as mentioned earlier, 
the Bianchi identities (4) guarantee that the evolution equations will preserve 
the constraints on future slices of the evolution. 

As we will see below, the Hamiltonian constraint (14) most naturally con- 
strains the 3-metric 7ij, while the momentum constraints (15) naturally con- 
strain the extrinsic curvature Kij . Taking the constraints into consideration, it 
seems that the 3-metric has five degrees of freedom remaining, while the extrin- 
sic curvature has three. But we know that the gravitational field in Einstein's 
theory has two dynamical degrees of freedom, so we expect that both and 
Kij should each have only two free components. The answer to this problem is, 
once again, the coordinate invariance of Einstein's theory. This seems strange 
at first, because we have already used the coordinate invariance of the theory 
to narrow our scope from the ten components of the 4-metric to the six compo- 
nents of the 3-metric. However, the lapse and shift do not completely specify the 
coordinate gauge. Rather, they specify how an initial choice of gauge will evolve 
with the foliation. The metric on a given hypersurface retains full 3-dimensional 
coordinate invariance, reducing the number of freely specifiable components to 
two [106; 107]. There also remains one degree of gauge freedom associated 
with the time coordinate which must be fixed. Each hypersurface represents a 
t = const, slice of the spacetime, so how the initial hypersurface is embedded in 
the full 4-dimensional manifold represents our temporal gauge choice. There is 
no unique way to specify this choice, but it is often convenient to let the trace 
of the extrinsic curvature K represent this temporal gauge choice [107]. Thus, 
we find that we are allowed to choose freely five components of the 3-metric and 

^The stress-energy tensor is decomposed as Tf^v = S^v + nnjv + Uvj^ H- n^rivp or, equiv- 
alents, Sij s T,j, ji = —n^^Ti^, and p = n^^n^Tfiv Note that the matter terms in Eqs. (12), 
(14), and (15) are defined with respect to a normal observer. This is in contrast to the usual 
definition of p with respect to the rest frame of the fluid in hydrodynamics. 
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three components of the extrinsic curvature. However, only two of the compo- 
nents for each field represent dynamical degrees of freedom, the remainder are 
gauge degrees of freedom. 

The four constraint equations, (14) and (15), represent conditions which 
the 3-metric and extrinsic curvature must satisfy. But, they do not specify 
which components (or combination of components) are constrained and which 
arc freely specifiable. In the weak field limit where Einstein's equations can be 
linearized, there are clear ways to determine which components are dynamic, 
which are constrained, and which are gauge. However, in the full nonlinear 
theory, there is no unique decomposition. In this case, one must choose a method 
for decomposing the constraint equations. The goal is to transform the equations 
into standard elliptic forms which can be solved given appropriate boundary 
conditions [80, 82, 110, 100]. Each different decomposition yields a unique set 
of elliptic equations to be solved and a unique set of freely specifiable parameters 
which must be fixed somehow. Seemingly similar sets of assumptions applied 
to different decompositions can lead to physically different initial conditions. 

2.2 York— Lichnerowicz Conformal Decompositions 

For general initial-data configurations, the most widely used class of constraint 
decompositions are the York-Lichnerowicz conformal decompositions. At their 
heart are a conformal decomposition of the metric and certain components of 
the extrinsic curvature, together with a transverse-traceless decomposition of 
the extrinsic curvature. 

First, the metric is decomposed into a conformal factor tp multiplying an 
auxiliary 3-metric [68, 106, 107]: 

lii=ip*lir (16) 

The auxiliary 3-metric jij is often called the conformal or background 3-metric, 
and it carries five degrees of freedom. Its natural definition is given by 

7ij=7"^ij, (17) 

leaving 7=1, but we are free to choose any normalization for 7. Using (16), 
we can rewrite the Hamiltonian constraint (14) as 

- - li'^K'^ + y^KijK'^ = -2ttG^^p, (18) 

where = V*Vi is the scalar Laplace operator, and and R are the covariant 
derivative and Ricci scalar associated with 7^^-. Equation (18) is a quasilinear 
elliptic equation for the conformal factor ip, and we see that the Hamiltonian 
constraint naturally constrains the 3-metric. 

The conformal decomposition of the Hamiltonian constraint was proposed 
by Lichnerowicz. But, the key to the full decomposition is the treatment of the 
extrinsic curvature introduced by York [108, 109]. This begins by splitting the 
extrinsic curvature into its trace and tracefree parts, 

Kij = Aij + l^ijK. (19) 
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The decomposition proceeds by using the fact that we can covariantly spUt any 
symmetric tracefree tensor as foUows: 



S'^ = {LXy^ +T*^. (20) 

Here, T'-^ is a symmetric, transverse-traceless tensor (i.e., VjT^^ = and T? = 0) 
and 

{hXy^ = V'X^ + V^X' - §7*^ V^X^. (21) 



After separating out the transverse-traceless portion of S^^ , what remains, {hxy^ , 
is referred to as its "longitudinal" part. We now want to apply this transverse- 
traceless decomposition to the tracefree part of the extrinsic curvature Aij. 
However, the conformal decomposition of the metric leaves us with at least two 

ways to proceed. 

The goal of the decomposition is to produce a coupled set of elliptic equa- 
tions to be solved with some prescribed boundary conditions. We have already 
reduced the Hamiltonian constraint to an elliptic equation being solved on a 
background space in terms of differential operators that are compatible with the 
conformal 3-metric. In the end, we want to reduce the momentum constraints 
to a set of elliptic equations based on differential operators that arc compatible 
with the same conformal 3-metric. But, the longitudinal operator (21) can be 
defined with respect to any metric space. In particular, it is natural to consider 
decompositions with respect to both the physical and conformal 3-metrics. 

2.2.1 Conformal Transverse- Traceless Decomposition 

Let us first consider decomposing A^^ with respect to the conformal 3-metric [110, 
115]. As we will see, when certain assumptions are made, this decomposition 
has the advantage of producing a simpler set of elliptic equations that must be 
solved. The first step is to define the conformal tracefree extrinsic curvature 
A'^ by 

A'^=ip-^°A'^ or Ay=V"^i,j. (22) 

Next, the transverse-traceless decomposition is applied to the conformal extrin- 
sic curvature, 

A'^ = {Lxy^ + Q'^ . (23) 

Note that the longitudinal operator L and the symmetric, transverse-tracefree 
tensor Q^^ are both defined with respect to covariant derivatives compatible 
with 7jj. 

Applying equations (16), (19), (21), (22), and (23) to the momentum con- 
straints (15), we find that they simplify to 

AlX' = lijj^V'K + SnGip^^f, (24) 

where 

AlX' = Vj(LXy^ = V^X' + \V\VjX^) + RiX\ (25) 
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and we have used the fact that 



Vj.S'^ =^-i"Vj(V''"5'^) (26) 

for any symmetric tracefree tensor S"^^ . 

In deriving equation (24), we have also used the fact tliat Q'^ is transverse 
(i.e. VjQ^^ = 0). However, in general, we will not know if a given symmet- 
ric tracefree tensor, say M*-', is transverse. By using (20) we can obtain its 
transverse-traceless part Q'-' via 

Qy=My_(Ly)y, (27) 

and using the fact that if Q*-' is transverse, we find 

VjQ'^ = = VjM'^ - Aj^YK (28) 

Thus, Eqs. (27) and (28) give us a general way of constructing the required 
symmetric transverse-traceless tensor from a general symmetric traceless tensor. 
Using the linearity of L, we can rewrite (23) as 

A'^ = {Lvy^ + (29) 

where 

V' = X'- Y\ (30) 
Similarly, using the linearity of Al, we can rewrite (24) as 

AlV* = Itp^V'K - VjM'^ + 87rGV'^°i'. (31) 

By solving directly for , wc can combine the steps of decomposing M'^ with 
that of solving the momentum constraints. 

After applying (19) and (22) to the Hamiltonian constraint (18), we obtain 
the following full decomposition, which I will list together here for convenience: 

la = tp%j 

A'^ = {Lvy^ + M'^ (32) 

V^i^-li^R-^i^^K^ + ^^-^AijA'^ = -27rGVV 

In the decomposition given by (32), we arc free to specify a symmetric tensor 
jij as the conformal 3- metric, a symmetric tracefree tensor M'-^, and a scalar 
function K. Then, with given matter energy and momentum densities, p and j% 
and appropriate boundary conditions, the coupled set of constraint equations 
for ip and are solved. Finally, given the solutions, we can construct the 
physical initial data, 'fij and K^^ . 
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The decomposition outlined above has the interesting property that if we 
choose K to be constant and if the momentum density vanishes'*, then the 
momentum constraint equations fully decouple from the Hamiltonian constraint. 
As we will see later, this simplification has proven to be useful. 

2.2.2 Physical Transverse- Traceless Decomposition 

Alternatively, we can decompose A'-' with respect to the physical 3-metric [81, 
82, 83]. We decompose the extrinsic curvature as 

A'^ = {Lwy^ + Q'^ . (33) 

In this case, the longitudinal operator L and the symmetric transverse-tracefree 
tensor Q*-' are both defined with respect to covariant derivatives compatible 
with 7y . 

Applying equations (16), (19), (21), (33), and (26) to the momentum con- 
straint (15), we find that it simplifies to 

AlW* + 6{Lwy^^j In V = IV'K + SnGij^f, (34) 

where we have used the fact that 

{Lwy^ = i}-'^{lwy^ . (35) 

As in the previous section, we will obtain the symmetric transvcrse-traceless 
tensor Q*-' from a general symmetric tracefree tensor M*-' by using (20). In this 
case, we take 

Q'^ = V-^°M^^' - {Lzy^, (36) 
and use the fact that Q*-' is transverse, to obtain 

Aj^Z' + 6{Lzy^Vj lnil) = V^VjM*^'. (37) 

Again, we can define 

V' = W'-Z\ (38) 

and use the linearity of L and Al to combine the process of obtaining the 
transverse-traceless part of M'-' and solving the momentum constraints. We 
obtain the following full decomposition, which I will list together here for con- 
venience^: 

A'^ = {Lvy^ +ij-^M'^ (39) 
AlF^ + 6{lvy^Wj InV' = fV^if - i^-^VjM'^ + STrGip^f 

*Or, if we solve the momentum constraint in terms of a ba<;kground momentum density 
s tp^'^j^, the momentum constraint decouples with nonvanishing j'. 

^Note that in (39) we have not used the usual conformal scaling of the extrinsic curvature 
given in equation (22). 
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In the decomposition given by (39) , we are again free to specify a symmetric 
tensor 7^- as the conformal 3-metric, a symmetric tracefree tensor M^^ , and a 
scalar function K. Then, with given matter energy and momentum densities, 
p and 7*, and appropriate boundary conditions, the coupled set of constraint 
equations for ^ and are solved. Finally, given the solutions, we can construct 
the physical initial data, jij and K'^^ . 

Notice that, while very similar to the decomposition from §2.2.1, the sets of 
equations are distinctly different. In general, if we make the same choices for 
the freely specifiable data in both decompositions (i.e., we choose -jij, M'-', and 
K the same), wc will produce two different sets of initial data. Both will be 
equally valid solutions of the constraint equations, but they will have distinct 
physical properties. 

There is at least one exception to this. Assume we have a valid set of initial 
data and Kij, which satisfies the constraint equations (14) and (15). For any 
everywhere-positive function ^, we define our freely specifiable data as follows: 



Then the solution to both sets of equations, assuming we use correct boundary 
conditions, will he ip = and V"^ = 0, which yields the original data as the 
solution for each decomposition. 

2.3 Thin- Sandwich Decomposition 

The two general initial- value decompositions outlined in §2.2.1 and §2.2.2 re- 
quire identical freely specified data (7^-, M'-^ , and K), yet they usually produce 
different physical initial data. One shortcoming of these approaches is that they 
provide no direct insight into how to choose the freely specifiable data. All of the 
data arc determined by schemes that involve only a single spacclike hypersurfacc. 
The resulting constraint equations are independent of the kincmatical variables 
a and that govern how the coordinates move through spacetime, and thus 
there is no connection to dynamics. York's thin-sandwich decomposition [114] 
takes a different approach by considering the evolution of the metric between 
two neighboring hypersurfaces (the thin sandwich). This decomposition is very 
similar to an approach originated by Wilson [103, 46], but is somewhat more 
general. Perhaps the most attractive feature of this decomposition is the insight 
it yields into the choice of the freely specifiable data. 

The decomposition begins with the standard conformal decomposition of the 
3-metric (16). However, we next make use of the evolution equation for the met- 
ric (13) in order to connect the 3-metrics on the two neighboring hypersurfaces. 
Label the two slices by t and t' , with t' = t + 5t, then 7^'^- = 7^- + (dtjij) St. We 
would like to specify how the 3-metric evolves, but we do not have full freedom 
to do this. We know we can freely specify only the conformal 3-metric, and 




(40) 



K 
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similarly, we arc free to specify only the evolution of the conformal 3-metric. 
We make the following definitions: 



My =7354(7 ^Jij), (41) 

Uij = d^ij, (42) 

and 

f'u,j = 0. (43) 

The latter definition is made for convenience, so that we can treat tp, ^ij, and Uij 
as regular scalars and tensors instead of as scalar- and tensor-densities within 
this thin-sandwich formalism. 

The conformal scaling of Uij follows directly from (16), (41), (42), (43), and 
the identity that, for any small perturbation, 51n7 = Y^S^ij. The result is 

Uij = ip'^Uij, (44) 

which relies on the useful intermediate result that^ 

da = 0. (45) 

Equation (41) represents the tracefree part of the evolution of the 3-metric, 
so (13) becomes 

u'^ = -2aA'^ + {Lf3Y^ . (46) 
Using the conformal scalings (22), (35), and (44), we obtain 

A'^ = ^ ((£/?)'•'■ - ff^^ . (47) 

York has pointed out that it is natural to use the following conformal rescal- 
ing of the lapse: 

a = Tp^a. (48) 

This rcscaling follows naturally from the "slicing function" that replaces the 
usual lapse (a = ^/7(5) which has been critical in solving several problems [4]. 
It also results in the natural conformal scaling (22) postulated for the tracefree 
part of the extrinsic curvature. Substituting (48) into (47) yields what is taken 
as the definition of the tracefree part of the conformal extrinsic curvature, 

A'^ = ^ ({Ipy^ - u'A . (49) 

Because the tracefree extrinsic curvature satisfies the normal conformal scal- 
ing, the Hamiltonian constraint will take on the same form as in (32). However, 

^Equations (43) and (45) take the place of defining 7=1 and the subsequent necessity of 
treating ip and 7ij as densities. 
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the momentum constraint will have a very different form. Combining equa- 
tions (16), (19), (21), (22), and (49) with the momentum constraint (15), we 
find that it simplifies to 

- {Ipy^^j lna = -^aip^V'K + 5V {lu'^) + 167rGai/'i°/. (50) 

Let us, for convenience, group together all the equations that constitute the 
thin-sandwich decomposition: 



^ ({Ipr V') (51) 



2a 

- {tfiy^Vj In 5 + ^aip^V'K = 5V (iu*^) + IQirGaip^^f 
V'^ip - li}R - ^ip^K^ + ^^-'^AijA'^ = -27rGV'^p 

In this decomposition (51), we are free to specify a symmetric tensor as the 
conformal 3-metric, a symmetric tracefree tensor ii^^ , a scalar function K, and 
the scalar function a. Solving this set of equations with appropriate boundary 
conditions yields initial data 7^ and Kij on a single hypersurface. However, 
we also know the following. // we chose to use the shift vector obtained from 
solving (50) and the lapse from (48) via our choice of a and our solution to the 
Hamiltonian constraint, then the rate of change of the physical 3-metric is given 

by 

dan = + flii (Vfc/3'= - aif ) (52) 
= [uij + 1% (Vfc/3'= + 6/3'= Vfe In V - V^ai^)' • 

This direct information about the consequences of our choices for the freely 
specifiable data is something not present in the previous decompositions. As we 
will see later, this framework has been used to construct initial data that are in 
quasiequilibrium. 



2.4 Stationary Solutions 

When there is sufficient symmetry present, it is possible to construct initial data 
that are in true equilibrium. These solutions possess at least two Killing vectors, 
one that is timelike at large distances and one that is spatial, representing an 
azimuthal symmetry. When these symmetries are present, solving for the initial 
data produces a global solution of Einstein's equations and the solution is said 
to be stationary. The familiar Kerr-Neumann solution for rotating black holes 
is an example of a stationary solution in vacuum. Stationary configurations 
supported by matter are also possible, but the matter sources must also satisfy 
the Killing symmetries, in which case the matter is said to be in hydrostatic 
equilibrium [8]. 
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The basic approach for finding stationary solutions begins by simphfying the 
metric to take into account the symmetries. Many different forms have been 
used for the metric (cf. Refs. [8, 9, 22, 34, 62, 21]). I will use a decomposition 
that makes comparison with the previous decompositions straightforward. First, 
define the interval as 

ds2 = + [A2(dr2 + r'^Ae'^) + sin^ ^(d<^ + fi'^dtf] . (53) 

This form of the metric can describe any stationary spacetime. Notice that the 
lapse is related to the conformal factor by 

a = (54) 

and that the shift vector has only one component 

/3' = (0,0,/?^). (55) 

I have used the usual conformal decomposition of the 3-metric (16) and have 
written the conformal 3-metric with two parameters as 

/ ^2 \ 

= A^r^ . (56) 

V SVsin^^ / 

The four functions if), (3'^, A, and B are functions of r and 6 only. 

The equations necessary to solve for these four functions are derived from 
the constraint equations (14) and (15), and the evolution equations (12) and 
(13). For the evolution equations, we use the fact that dtjij = and dtKij = 0. 
The metric evolution equation (13) defines the extrinsic curvature in terms of 
derivatives of the shift 

^..^^(Vi/?,+V,A). (57) 

With the given metric and shift, we find that K = and the divergence of the 
shift also vanishes. This means we can write the tracefree part of the extrinsic 
curvature as 

^ii = ^-ioj[^j = ^^p~^{li3Yi. (58) 

We find that the Hamiltonian and momentum constraints take on the forms 
given by the thin-sandwich decomposition (51) with u'-' = K = and a = 
Only one of the momentum constraint equations is non-trivial, and we find that 
the constraints yield elliptic equations for ip and f3'^. What remains unspecified 
as yet are A and B (i.e., the conformal 3-metric). 

The conformal 3-metric is determined by the evolution equations for the 
traceless part of the extrinsic curvature. Of these five equations, one can be 
written as an elliptic equation for B, and two yield complementary equations 
that can each be solved by quadrature for A. The remaining equations are 
redundant as a result of the Bianchi identities. 

Of course, the clean separation of the equations I have suggested above 
is an illusion. All four equations must be solved simultaneously, and clever 
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combinations of the four metric quantities can greatly simplify the task of solving 
the system of equations. This accounts for the numerous different systems used 
for solving for stationary solutions. 
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3 Black Hole Initial Data 



In this section, we will look at Cauchy initial data that represent one or more 
black holes in an asymptotically flat spacctime. The majority of these will be 
either vacuum solutions or solutions of the Einstein-Maxwell equations. With 
no matter to support the gravitational field, we find that we must usually^ use 
a spacctime with a non-trivial topology. It is certainly possible to construct 
black- hole solutions supported by matter [89, 90, 5], but it is often desirable to 
avoid the complications of matter sources. 

This raises a point about solutions of Einstein's equations which we have 
not yet mentioned. When constructing solutions of Einstein's initial- value equa- 
tions, we are free to specify the topology of the initial-data hypersurface. Ein- 
stein's equations of general relativity place no constraints on the topology of 
the spacctime they describe or of spacelike hypersurfaces that foliate it. For 
astrophysical black holes (i.e., black holes in an asymptotically flat spacetime), 
the freedom in the choice of the topology has relatively minor consequences. 
The primary effects of different topology choices are hidden within the black 
hole's event horizon. 

In the sections below, we will explore many of the existing black-hole solu- 
tions and the schemes for generating them. 

3.1 Classic Solutions 
3.1.1 Schwarzschild 

Of course, the simplest black-hole solution is the Schwarzsc;liild solution. It 
represents a static spacctime containing a single black hole that connects two 
causally disconnected, asymptotically flat universes. There are actually many 
different coordinate representations of the Schwarzschild solutions. The simplest 
representations are time-symmetric {Kij = 0), and so exist on a "maximally 
embedded" spacelike hypersurface (K = 0). These choices fix the foliation {E}. 
Spherical symmetry fixes two of the three spatial gauge choices. If we choose 
an "areal-radial coordinate"^, then the interval is written as 



= _ (1 _ 2M) + (1 _ 2M)-i dr^ + r^{d0^ + sin^ Odc^^). (59) 



If we choose an isotropic radial coordinate, then the interval is written as 



In both (59) and (60), M represents the mass of the black hole as measured 
at spacelike infinity. Both of these solutions exist on the same foliation of 
t = const, slices. But, notice that the 3-geometry of the slice associated with 
(60) is conformally fiat, while the 3-geometry associated with (59) is not. 

'^A black hole can be supported by a compact gravitational wave [2]. 
^Spheres at constant areal-radial coordinate r have proper area 47rr^. 




(60) 
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The solution given in (60) is easily generated by any of the methods in §2.2 or 
§2.3. By choosing a time-symmetric initial-data hypersurface, we immediately 
get Kij = 0, which eliminates the need to solve the momentum constraints. If 

we choose the conformal 3-geometry to be given by a flat metric (in spherical 
coordinates in this case), then the vacuum Hamiltonian constraint (18) becomes 

V^V = 0, (61) 

where is the flat-space Laplace operator. For the solution tp to yield an 
asymptotically flat physical 3-metric, we have the boundary condition that 
■^(f — > oo) = 1. The simplest solution of this equation is 

M 

^ = 1 + -, (62) 

where we have chosen the remaining integration constant to give a mass at 
infinity of M. 

We now have full Cauchy initial data representing a single black hole. If we 

want to generate a full solution of Einstein's equations, we must choose a lapse 
and a shift vector and integrate the evolution equations (12) and (13). In this 
reasonable approach for specifying the lapse is to demand that the time 
derivative of K vanish. For the case of K = 0, this yields the so-called maximal 
slicing equation which, for the current situation, takes the form 

V^imJj) = 0. (63) 

If we choose boundary conditions so that the lapse is frozen on the event horizon 
(a(f = M/2) = 0) and goes to one at infinity, we find that the solution is 

1 _ M 

"=7-f- (64) 

If we now choose /?* = 0, we find that the left-hand sides of the evolution 
equations (12) and (13) vanish identically, and we have found the static solution 
of Einstein's equations given in (60). We can, of course, recover the usual 
Schwarzschild coordinate solution (59) by using the purely spatial coordinate 
transformation r = r(l -|- ^)^. 

It is interesting to examine the differences in these two representations of 
the Schwarzschild solution. The isotropic radial coordinate representation is 
well behaved everywhere except, it seems, at r = 0. However, even here, the 
solution is well behaved. The 3-geometry is invariant under the coordinate 
transformation 

'"^^ ' (65) 

The event horizon at f = is a fixed-point set of the isometry condition 
(65) which identifies points in two causally disconnected, asymptotically flat 
universes. We see that f = is simply an image of infinity in the other uni- 
verse [31]. 
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Given our choice for the lapse (64), which is frozen on the event horizon, 
we find that the solution can cover only the exterior of the black hole. To 
cover any of the interior with the lapse pinned to zero at the horizon would 
require wc use a slice that is not spacclike everywhere. This is exactly what 
happens when the usual Schwarzschild areal-radial coordinate is used. At the 
event horizon, r = 2M, there is a coordinate singularity, and inside this radius 
the t = const, hypcrsurfacc is no longer spacclike. It is impossible to perform a 
Cauchy evolution interior to the event horizon using the areal-radial coordinate 
and the given time slicing. 

We find that a Cauchy evolution, using the usual Schwarzschild time slicing 
that is frozen at the horizon, is capable of evolving only the region exterior to 
the black hole's event horizon. Portions of the interior of the black hole can be 
covered by an evolution that begins with data on a standard Schwarzschild time 
slice, but the result is not a time-independent solution. As we will see later, 
there are other slicings of the Schwarzschild spacetime that cover the interior of 
the black hole and yield time-independent solutions. 



3.1.2 Time- Symmetric Multi-Hole Solutions 

As we saw in §3.1.1, the simplest approach for generating initial data is to assume 
time symmetry and let the conformal 3-geometry be flat. Wc could have gen- 
erated the exterior solution for Schwarzschild with the areal-radial coordinate 
by a clever choice for the conformal 3-geometry. The Hamiltonian constraint is 
still linear, even if jij is not flat. But there is no obvious motivation for the 
correct that would yield the desired solution. 

One approach for generating a time-symmetric multi-hole solution is straight- 
forward. Brill Lindquist initial data [31, 69] again assume a flat conformal 
3-geometry, and the only non-trivial constraint equation is the Hamiltonian 
constraint, which again takes the form given in (61). But this time, we use 
the linearity of the Hamiltonian constraint and choose the solution to be a su- 
perposition of solutions with the form of (62). More precisely, we choose the 
solution 

N 

Here, |x — C(j| is a coordinate distance from the point in the Euclidean 
conformal space, and the /x^ are constants related to the masses of the black 
holes. Assuming the points are sufficiently far apart, this solution of the 
initial-value equations represents black holes momentarily at rest in "our" 
asymptotically flat universe. As was the case for a single black hole, each sin- 
gular point in the solution, x = C^, represents infinity in a different, causally 
disconnected universe. In fact, each black hole connects "our" imiversc to a 
different universe, so that there are A^ -|- 1 asymptotically flat hypersurfaces 
connected together at the throats of N black holes. While we started with a 
3-dimensional Euclidean manifold, the requirement that we delete the singular 
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points, C(j, results in a manifold that is not simply connected. This solution is 
often referred to as having a topology with + 1 "sheets" . 

Brill-Lindquist initial data are very similar to the Schwarzschild initial data 
in isotropic coordinates except for one major difference: the solution docs not 
represent two identical universes that have been joined together. The coordi- 
nate transformation (65) can still be used to show that each pole in the solution 
corresponds to infinity on an asymptotically flat hypersurface. But, the solution 
has N + 1 different asymptotically flat universes connected together, not two, 
and "our" universe containing N black holes cannot be isometric to any of the 
other universes which contain just one. Interestingly, Misncr [78] found that it 
is possible to construct a solution of the vacuum, time-symmetric Hamiltonian 
constraint (61) that has two isometric asymptotically flat hypersurfaces con- 
nected by N black holes. The case of two black holes that satisfies this isomctry 
condition (often called inversion symmetry) is usually referred to as "Misner 
initial data". It has an analytic representation in terms of an infinite series 
expansion. The construction is tedious, and there are several representations of 
the solution [92, 38]^. 

Like the Brill-Lindquist data, the non-simply connected topology of the full 
manifold is represented on a Euclidean 3-manifold by the presence of singiilar 
points that must be removed. Brill-Lindquist data have N singular points, each 
representing an image of infinity as seen through the throat of the black hole 
connecting our universe to that hole's other universe. Misner's data contain 
an infinite number of singular points for each black hole, each representing 
an image of one of two asymptotic infinities. The result is seen to represent 
two identical asymptotically flat universes joined by N black holes. The two 
universes join together at the throats of the N black holes, with each throat 
being a coordinate 2-sphere in the conformal space. Data at any point in one 
universe are related to data at the corresponding point in the alternate universe 
via the same isometry condition (65) found in the Schwarzschild case. And, 
as in the Schwarzschild case, the 2-sphere throat of each black hole forms a 
flxed-point set of the isometry condition. 

Both the Brill-Lindquist data and Misner's data represent N black holes 
at a moment of time-symmetry (i.e., all the black holes are momentarily at 
rest). Both are conformally flat and the difference in the topology of the two 
solutions is hidden from an observer outside the black holes. Yet solutions 
where the holes are chosen to have the same size and separation yield similar, 
but physically distinct, solutions [36, 3]. 

3.2 General Multi-Hole Solutions 

Time-symmetric black-hole solutions of the constraint equations such as those 
described in §3.1 are useful as test cascis bc;causc! they have analytic repre- 
sentations. However, they have very little physical relevance. General time- 
asymmetric solutions are needed to represent black holes that are moving and 

^In Ref. [38], Eq. (B7), the line "1 for n = 1" should read "a"! for n = 1". 
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spinning"'^'' . A few approaches for generating general multi-hole solutions have 
been explored and, below, we will look at the two approaches which are direct 
generalizations of the Misner and Brill-Lindquist data. Generalizing these two 
approaches was the most natural first step toward constructing general multi- 
hole solutions. 

The first approach to be developed generalized the Misner approach [78] . It 

was attractive because an isomctry condition relating the two asymptotically 
flat universes provides two useful things. First, because the two universes are 
identical, finding a solution in one universe means that you have the full solution. 
Second, because the throats arc fixed-point sets of the isometry, we can construct 
boundary conditions on any quantity there. This allows us to excise the region 
interior to the spherical throats from one of the Euclidean background spaces 
and solve for the initial data in the remaining volume. 

The generalization of the Misner approach seemed preferable to trying to 
solve the constraints on N +1 Euclidean manifolds stitched together smoothly 
at the throats of N black holes. However, Brandt and Briigmann [27] realized 
it was possible to factor out analytically the behavior of the singular points 
in the Euclidean manifold of the N + 1 sheeted approach. Referred to as the 
"puncture" method, this approach allows us to rewrite the constraint equations 
for functions on an -|- 1 sheeted manifold as constraint equations for new 
functions on a simple Euclidean manifold. 

Another approach tried, which we will not discuss in detail, avoided the issue 
of the topology of the initial-data slice entirely. Developed by Thornburg [99], 
this approach was based on the idea that only the domain exterior to the ap- 
parent horizon of a black hole is relevant. The equation describing the location 
of an apparent horizon can be rewritten in a form that can be used as a bound- 
ary condition for the conformal factor in the Hamiltonian constraint equation. 
Thus, given a compatible solution to the momentum constraints, this bound- 
ary condition can be used to construct a solution of the Hamiltonian constraint 
in the domain exterior to the apparent horizons of any black holes, with no 
reference at all to the topology of the full manifold. 

3.2.1 Bowen-York Data 

The generalization of the Misner approach was developed by York and his col- 
laborators [26, 23, 24, 66, 112, 65, 25, 38]. The approach is often called the 
"conformal-imaging" method, and the data are usually referred to as "Bowen- 
York" data. This approach begins with a set of simplifying assumptions that is 
common to all three of the approaches described above. These assumptions are 
that 

K = 0, m,aximal slicing 

lij — fij I conformal flatness (67) 

V'loo = 1. asymptotic flatness 

1" Time-asymmetric solutions are also needed to represent time-independent solutions that 
cover the interior of a black hole. 
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Here, fij represents a flat metric in any suitable coordinate system. The assump- 
tion of conformal flatness means that the differential operators in the constraints 
are the familiar flat-space operators. More importantly, if we use the conformal 
transverse-traceless decomposition (32), we find that in vacuum the momentum 
constraints completely decouple from the Hamiltonian constraint. 

The importance of this last property stems from the fact that York and 
Bowen were able to find analytic solutions of this version of the momentum 
constraints, solutions that represent a black hole with both linear momentum 
and spin [26, 23, 24]. If we choose M*^ = 0, then the momentum constraints 
(31) become 

+ iV* (VjV^) = 0. (68) 

A solution of this equation is 

V' = [7P' + n'rijP'] + ^t'^^UjSk. (69) 

Here, and S"^ are vector parameters, r is a coordinate radius, and is the 
outward-pointing unit normal of a sphere in the flat conformal space (n* = ^). 
e*-^*^ is the 3-dimensional Levi-Civita tensor. 

This solution of the momentum constraints yields the tracefree part of the 
extrinsic curvature 

^io = -^[Pinj+Pjni-{fij-ninj)P''nk] (70) 

Remarkably, using this solution (70) and the assumptions in (67), we can de- 
termine the physical values for the linear and angular momentum of any initial 

data we can construct. The momentum contained in an asymptotically flat 
initial-data hypersurface can be calculated from the integral [111] 

^%) = ^ £ {Ki - SiK) q.^d'S,, (71) 

where ^^^^ is a Killing vector of the 3- metric 7^^^. Since we are not likely 
to have true Killing vectors, we make use of the asymptotic translational and 
rotational Killing vectors of the flat conformal space. We find from (71), (70), 
and (67) that the physical linear momentum of the initial-data hypersurface is 
P' and the physical angular momentum of the slice is 5'. Furthermore, because 
the momentum constraints are linear, we can add any number of solutions of 
the form of Eq. (70) to represent a collection of linear and angular momentum 
sources. The total physical linear momentum of the initial-data slice will simply 

^^If is a translational Killing vector, then (71) yields the linear momentum in the 
direction of that Killing vector. If ^^^^j is a rotational Killing vector, then (71) yields the 
corresponding angular momentum. 
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be the vector sum of the individual hnear momenta. The total physical angular 
momentum cannot be obtained by simply summing the individual spins because 
this neglects the orbital angular momentum of the various sources. However, 
the total angular momentum can still be computed without having to solve the 
Hamiltonian constraint [113]. 

The Bowen-York solution for the extrinsic curvature is the starting point for 
all the general multi-hole initial-data sets we have discussed in §3.2. However, 
this solution is not inversion symmetric. That is, it does not satisfy the isometry 
condition that any field must satisfy to exist on a two-sheeted manifold like that 
of Misncr's solution. Fortimatcly, there is a method of images, similar to that 
used in electrostatics but applicable to tensors, that can be used to make any 
tensor inversion symmetric [78, 26, 66, 65, 112]. 

For the conformal extrinsic curvature of a single black hole, there are two 
inversion-symmetric solutions [26] 

3 

T ^ [Ptrij + PjUr + iftj - 5ninj)P''nk] (72) 
3 

+ ^ [ckieS^n'^nj + CkjeS^n'^ni] . 

Here, a is the radius of the coordinate 2-sphere that is the throat of the black 
hole. Of course, this coordinate 2-sphere is the fixed-point set of the isometry 
and is the surface on which we can impose boundary conditions. Notice that 
this radius enters the solutions only when we make it inversion symmetric. 

When the extrinsic curvature represents more than one black hole, the pro- 
cess for making the solution inversion symmetric is rather complex and results in 
an infinite-series solution. However, in most cases of interest, the solution con- 
verges rapidly and it is straightforward to evaluate the solution numerically [38] . 

Given an inversion-symmetric conformal extrinsic curvature, it is possible to 
find an inversion-symmetric solution of the Hamiltonian constraint [26]. Given 
our assumptions (67), the Hamiltonian constraint becomes 

VV + li'-^AijA'^ = 0. (73) 

The isometry condition imposes a condition on the conformal factor at the 
throat of each hole. This condition takes the form [26] 



2ra 



(74) 



where is the outward-pointing unit-normal vector to the a^^ throat and is 
the coordinate radius of that throat. This condition can be used as a boundary 
condition when solving (73) in the region exterior to the throats. 

In addition to boundary conditions on the throats, a boundary condition 
on the outer boundary of the domain is needed before the quasilinear elliptic 
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equation in (73) can be solved as a well-posed boundary-value problem. This 
final boundary condition comes from the fact that we want an asymptotically 
flat solution. This implies that the solution behaves as 

V' = l + |+0(r-^), (75) 

where E is the total ADM energy content of the initial-data hypersurface. Equa- 
tion (75) can be used to construct appropriate boundary conditions either at 
infinity or at a large, but finite, distance from the black holes [115]. 



3.2.2 Puncture Data 

The generalization of the Brill- Lindquist data developed by Brandt and Briig- 
mann [27] begins with the same set of assumptions (67) as the conformal-imaging 
approach outlined in §3.2.1. We immediately have (70) from the solution of the 
momentum constraints, and we must solve the Hamiltonian constraint, which 
again takes the form of Eq. (73). At this point, however, the method of solution 
differs from the conformal-imaging approach. 

Based on the time-symmetric solution, it is reasonable to assume that the 
conformal factor will take the form 

^ = ^+u, - = y , , . (76) 

X X ^,2|x-C,| ^ ' 

If u is sufficiently smooth, (76) implies that the manifold will have the topology 
of + 1 asymptotically flat regions just as in the Brill-Lindquist solution. In 
this case, asymptotic flatness requires that u = 1 + 0{r~^). 
Substituting (76) into the Hamiltonian constraint (73) yields 

V^u + r]{l + xu)-^ = 0, (77) 

where 

V=lx'AijA''- (78) 

Near each singular point, or "pimcture", we find that x « |x— C^l. From (70), 
we see that AijA^^ behaves no worse than |x — C<t|~®, so ij vanishes at the 
punctures at least as fast as |x — C^l- 

With this behavior, Brandt and Briigmann [27] have shown the existence 
and uniqueness of solutions of the modified Hamiltonian constraint (77). 
The resulting scheme for constructing multiple black hole initial data is very 
simple. The mass and position of each black hole are parameterized by /io- and 
Co-, respectively. Their linear momenta and spin are parameterized by and 
So- in the conformal extrinsic curvature (70) used for each hole. Finally, the 
solution for u is found on a simple Euclidean manifold, with no need for any 
inner boundaries to avoid singularities. This is a great simplification over the 
conformal-imaging approach, where proper handling of the inner boundary is 
the most difficult aspect of solving the Hamiltonian constraint numerically [41]. 
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3.2.3 Problems with These Data 



Both the conformal-imaging and puncture methods for generating muhiple black 
hole initial data allow for completely general configurations of the relative sizes 
of the black holes, as well as their linear and angular momenta. This does not 
mean that these schemes allow for the generation of all desired black- hole initial 
data. The two schemes rely on specific assumptions about the freely specifiable 
gravitational data. In particular, they assume K = 0, M'-' = 0, and, most 
importantly, that the 3-gcomctry is conformally flat. 

These choices for the freely specifiable data are not always commensurate 
with the desired physical solution. For example, if we choose to use either 
method to construct a single spinning black hole, we will not obtain the Kerr 
solution. The Kerr-Newman solution can be written in terms of a quasi-isotropic 
radial coordinate on a K = time slice [28]. Let r denote the usual Boyer- 
Lindquist radial coordinate and make the standard definitions 



+ cos^ e and A = - 2Mr + + 



(79) 



A quasi-isotropic radial coordinate f can be defined via 



1 + 



M + + Q2 
2f 



1 + 



2f 



(80) 



The interval then becomes 



e2^/3(d-2 ^ -2^^2) ^ -2 ^.^2 ^e-4^/3(d0 + p'^dtf] , (81) 



with 



a 



p2A 



+ - A 



= —a 



(r2 + a2)2 - Aa2 sin^ 9 
/93 ((r^ + a^)^ - Aa^ sin^ 61) : 

^2 

(r2+a2)2- Aa^sin^^ 



(82) 



Wc see immediately that the 3-geometry associated with a t — const, hy- 
persurface of (81) is not conformally flat. In fact, Garat and Price [55] have 
shown that in general there is no spatial slicing of the Kerr spacetime that is 
axisymmetric, conformally flat, and smoothly goes to the Schwarzschild solution 
as the spin parameter a ^ 0. 

Since the Kerr solution is stationary, the inescapable conclusion is that con- 
formally flat initial data for a single rotating black hole must also contain some 
nonvanishing dynamical component. When we evolve such data, the system 
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will emit gravitational radiation and eventually settle down to the Kerr geome- 
try [26, 29]. But, it cannot be the Kerr geometry initially, and it is unlikely that 
the spurious gravitational radiation content of the initial data has any desirable 
physical properties. Conformally flat initial data for spinning holes contain some 
amount of unphysical "junk" radiation^^. 

The choice of a conformally flat 3-geometry was originally made for con- 
venience. Combined with the choice of maximal slicing, these simplifying as- 
sumptions allowed for an analytic solution of the momentum constraints which 
vastly simplified the process of constructing black-hole initial data. Yet there 
has been much concern about the possible adverse physical effects that these 
choices (especially the choice of conformal flatness) will have in trying to study 
black- hole spacetimes [40, 70, 87, 60, 86]. While these conformally flat data sets 
may still be useful for tests of black- hole evolution codes, it is becoming widely 
accepted that the unphysical initial radiation will significantly contaminate any 
gravitational waveforms extracted from evolutions of these data. In short, these 
data are not astrophysically realistic. 

The various initial-data decompositions outlined in §2.2 and §2.3 are all ca- 
pable of producing completely general black-hole initial data sets. The only 
limitation of these schemes is our understanding of what choices to make for 
the freely specifiable data and the boundary conditions to apply when solving 
the sets of elliptic equations. All these choices will have a critical impact on the 
astrophysical significance of the data produced. It is also important to remem- 
ber that similar choices for the freely specifiable data will result in physically 
different solutions when applied to the different schemes. 

The first studies of black-hole initial data that arc not conformally flat were 
carried out by Abrahams et al. [1]. They looked at the superposition of a 
gravitational wave and a black hole. Using a form of the conformal metric that 
allows for so-called Brill waves, they constructed time-symmetric initial data 
that were not conformally flat and yet satisfied the isometry condition (65) used 
in the conformal-imaging method of §3.2.1. These data were further generalized 
by Brandt and Seidel [30] to include rotating black holes with a superimposed 
gravitational wave. In this case, the data are no longer time-symmetric yet they 
satisfy a generalized form of the isometry condition so that the solution is still 
represented on two isometric, asymptotically flat hypersurfaces. 

Matzner et al. [77] have begun to move beyond conformally flat initial data 
for binary black holes^^. Their proposal is to use boosted versions of the Kerr 
metric written in the Kerr-Schild form to represent each black hole. Thus, an 
isolated black hole will have no spurious radiation content in the initial data. To 
construct solutions with multiple black holes, they propose, essentially, to use a 
linear combination of the single-hole solutions. The resulting metric can be used 
as the conformal 3-metric, the trace of the resulting extrinsic curvature can be 
used for K, and the tracefree part of the resulting extrinsic curvature can be used 

^^A similar conclusion is reached for conformally flat data for a single black hole with linear 
momentum [111]. 

'^■^A related method has been proposed by Bishop et al. [15], but their approach is much 
different and outside the current scope of this review. 
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for NP^ . Their scheme uses York's conformal transverse-traceless decomposition 
outUned in §2.2.1, with the boundary conditions of ^ = 1 and = on the 
horizons of the black holes and conditions appropriate for asymptotic flatness 
at large distances from the holes. 

The approach outlined by Matzner should certainly yield "cleaner" data than 
the conformally flat data currently available. For the task of specifying data for 
astrophysical black-hole binaries in nearly circular orbits, it is still true that this 
new data will not contain the correct initial gravitational wave content. Because 
the black holes are in orbit, they must be producing a continuous wave-train 
of gravitational radiation. This radiation will not be included in the method 
proposed by Matzner et al. Also, it is clear that the boundary conditions being 
used do not correctly account for the tidal distortion of each black hole by its 
companion. When the black holes arc sufficiently far apart, the radiation from 
the orbital motion can be computed using post-Newtonian techniques. One 
possibility for producing astrophysically realistic, binary black-hole initial data 
is to use information from these post-Newtonian calculations to obtain better 
guesses for 'jij, M^^ , and K. 



3.3 Horizon-Penetrating Solutions 

We noted in §3.1.1 that the time-independent maximal slicing of Schwarzschild 
with isotropic coordinates covers only the exterior of the black hole. This is 
because the time independence of this gauge requires that the lapse vanish on 
the horizon. It is possible to evolve into the black hole's interior when starting 
from initial data constructed in this gauge, but it requires a choice for the lapse 
that yields a time-dependent solution [93]. The time dependence of such a 
solution is purely gauge, of course, since the spacetime is static. 

It is possible to cover all, or part, of the interior of a single black hole with 
a time-independent slicing. However, doing so seems to require that we give up 
the maximal-slicing condition. To cover the interior of the black hole, we need 
a slicing that passes smoothly through the event horizon. A convenient way 
to generate such solutions is to begin with the metric in standard ingoing-null 
coordinates. If we want to consider a rotating and charged black hole, then we 
use the Kerr-Newman geometry in Kerr coordinates: 

ds^ = -fl- ^Mr-Q^ \ ^ 2dydr - 2 ~ a sin^g d^d^ (83) 



where p and A are defined by Eq. (79) and V is the ingoing-null coordinate. 
This metric is regular at r = r± = M ± — — Q'^, where r+ and r_ are 

the locations of the event horizon and the Cauchy horizon, respectively. 

This metric can be put into a form suitable for producing time-independent 
Cauchy initial data by making coordinate transformations of the general form 

dt = dV + /(r)dr, d^ = d^ -h g{r)dr, (84) 
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where / and g arc siiitably chosen functions of the radial coordinate r. There are 
a few particularly significant solutions for the general Kerr-Newman geometry, 
and I will outline these below, listing the nonzero components of the metric in 
the standard 3+1 format. 

3.3.1 Kerr— Schild Coordinates 

A spherical coordinate version of the standard Kerr-Schild coordinate system 
is obtained from (83) by using the coordinate choice (cf. Refs. [73, 42]) 

t = V-r and = 0. (85) 

The nonzero components of the lapse, shift, and 3- metric are then given by; 

2Mr - 



a 



1 + 



13'' = a 



2 2Mr - Q2 



1 + 



2Mr - Q2 



1 + 



2Mr - 0^ 



asin^9, 



199 = P , 



r +a -\ 7i a sm 6 



sin 



(86) 

(87) 

(88) 

(89) 
(90) 
(91) 



Cartesian coordinate components can be obtained from these via the standard 
Kerr-Schild coordinate transformations [79] 

X + iy = {r + ia)e^'^ sin 9 and z = r cos 6. (92) 

This yields the implicit definition of r from 

-r\x^+y^ + z'-a')-a^z^ = 0, (93) 

with r > and r = on the disk described by 2; = and x"^ +y'^ < a^. 

3.3.2 Harmonic Coordinates 

Harmonic time slicing is integral to some hyperbolic formulations of general rel- 
ativity, and a time-independent harmonic slicing of the Kerr-Newman geometry 
does exist [17, 42]. The harmonic time slicing condition is Ot = 0, which can 
be written 

-^d^iV^g^^) = 0. (94) 
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This equation is satisfied by using the coordinate choice 



t^V -r + 2M\n 



2M 



r — r- 



(95) 



The nonzero components of the lapse, shift, and 3-metric are then given by: 

a'^ = 1 + 



2Mr - Q2 


fr + r+\ 




f 2M \ 




\r — r- J 




— T- J 



p'' = a 



Irr — 
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2M 

r ~ 
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■'r+ 



r — r- 



■r+ 



1 + 



2Mr-Q'^ fr + r+ 



r — r- 



100 = P 

2 2Mr - 2 . 2, 
a H ^ a sm 6 



r — r- 



sm 



(96) 

(97) 

(98) 

(99) 

(100) 
(101) 

(102) 



Cartesian coordinate components can be obtained from these via the standard 
Kerr-Schild coordinate transformations (92) and (93). However, for the har- 
monic slicing, the t = const, hypersurface is spacelike only outside the Cauchy 
horizon at r > r_ . 

Fully harmonic coordinates (□a;'' = 0) can be defined when Cartesian spa- 
tial coordinates are used by employing a variation of the standard Kerr-Schild 
coordinate transformations [42] 

x + iy = {r — m + ia)e^'^ sin 6 and z = [r — m) cos 6. (103) 

This yields the implicit definition of r from 

(r - mf - (r - mf{x'^ + y"^ + - a^) - c? = 0. (104) 

Fully harmonic coordinates are useful because applying a boost to a harmoni- 
cally sliced black hole yields a solution that satisfies (94) only if the black hole 
is written in fully harmonic coordinates. In this case, the boosted solution also 
satisfies the fully harmonic coordinate conditions. 
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3.3.3 Generalized Painleve— GuUstrand Coordinates 



The Painlcvc GuUstrand gauge choice for the Schwarzschild geometry has been 
rediscovered many times because of its simple form (cf. Refs. [85, 59, 88, 64, 67]). 
It is another time-independent solution, but the 3-geometry is completely flat 
(not simply conformally flat). The lapse is one in this gauge, and all of the 
information regarding the curvature of spacetime is contained in the shift. The 
Painleve-GuUstrand gauge also has an intuitive physical interpretation [74]. 
An observer starting at rest at infinity and freely falling will trace out a world 
line that is everywhere orthogonal to the t = const, hypersurfaces in Painleve- 
GuUstrand coordinates. 

A generalization of the Painleve-GuUstrand gauge derived by Doran [48] 
includes the Kerr spacetime, and the extension of this solution to the full Kerr- 
Newman spacetime is trivial. In the limit that a and Q vanish, this solution 
reduces to the Painleve-GuUstrand gauge. The coordinate transformation is 
written most easily as 



dt = dV 



dr 



1 



2Mr-Q^ 



d4> = d^- 



dr 



1 + 



l 2Mr-Q^ 



(105) 



The nonzero components of the lapse, shift, and 3-metric are then given by: 

= 1 (106) 



a 



I 2Mr - Q2 



7r 



l2Mr-Q^ 



osm 
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(107) 
(108) 

(109) 

(110) 
(111) 



Notice that the lapse remains one, but the 3-gcomctry is no longer flat when 
the black hole is spinning. Cartesian coordinate components can be obtained 
from these via the standard Kerr-Schild coordinate transformations (92) and 
(93). Like the Kerr-Schild time slicing, a t = const, slice of the generalized 
Painleve-GuUstrand gauge remains spacelike for all r > 0. 
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3.4 Quasicircular Binary Data 

One of the primary driving forces behind the development of black-hole initial 
data has been the two-body problem of general relativity: the inspiral and co- 
alescence of a pair of black holes. This problem is of fundamental importance. 
Not only is the relativistic two-body problem the most fundamental dynamical 
problem of general relativity, it is also considered one of the most likely candi- 
dates for observation with the upcoming generation of gravitational wave laser 
interferometers. Because of the circularizing effects of gravitational radiation 
damping, we expect the orbits of most tight binary systems to have small ec- 
centricities. It is therefore desirable to have a method that can discern which 
data sets, within the very large parameter space of binary black-hole initial-data 
sets, correspond to black-hole binaries in a nearly circular (quasicircular) orbit. 

Currently, only one approach has been developed for locating quasicircular 
orbits in a parameter space of binary black- hole initial data [39] . It is based on 
the fact that minimizing the energy of a binary system while keeping the orbital 
angular momentum fixed will yield a circular orbit in Newtonian gravity. The 
idea docs not hold strictly for general relativistic binaries since they emit gravi- 
tational radiation and cannot be in equilibrium. However, for orbits outside the 
innermost stable circular orbit, the gravitational radiation reaction time scale is 
much longer than the orbital period. Thus it is a good approximation to treat 
such binaries as an equilibrium system. Called an "effective potential method" , 
this approach was used originally to find the quasicircular orbits and innermost 
stable circular orbit (ISCO) for equal-sized nonspinning black holes [39]. In this 
work, the initial data for binary black holes were computed using the conformal- 
imaging approach outlined in §3.2.1. The approach was also applied to binary 
black-hole data computed using the puncture method [10], where similar results 
were found. Configurations containing a pair of equal-sized black holes with 
spin also have been examined [86] . In this case, the spins of the black holes are 
equal in magnitude, but arc aligned cither parallel to, or anti-parallel to, the 
direction of the orbital angular momentum. 

The approach defines an "eff'ective potential" based on the binding energy 
of the binary. The binding energy is defined as 

Eb = Eat^m - Ml - M2, (112) 

where Eaum is the total ADM energy of the system measured at infinity, and 

Ml and M2 are the masses of the individual black holes. Quasicircular orbital 
configurations are obtained by minimizing the effective potential (defined as 
the nondimensional binding energy Eb/jJ, (where ^ = MiM2/(Mi -|- M2)) as a 
function of separation, while keeping the the ratio of the masses of the black 
holes M1/M2, the spins of the black holes Si/Mf and S2/MI, and the total 
angular momentum J/M1M2 constant. 

This approach is limited primarily by the ambiguity in defining the individual 
masses of the black holes. Mi and M2. There is no rigorous definition for the 
mass of an individual hole in a binary configuration and some approximation 
must be made here. There is also no rigorous definition for the individual spins of 
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holes in a binary configuration. The problem of defining the individual masses 
becomes particularly pronounced when the holes are very close together (see 
Ref. [86]). The limiting choice of conformal flatness for the 3-geometry also has 
proven to be problematic. The effects of this choice have been clearly seen in the 
case of quasicircular orbits of spinning black holes [86] , but it is also believed to 
be a serious problem for any binary configuration because binary configurations 
are not conformally flat at the second post- Newtonian order [87] . 

To date, the results of the effective-potential method have not matched well 
to the best result from post-Newtonian approximations [47]^^. It will be inter- 
esting to see if the results from post-Newtonian approximations and numerical 
initial-data sets converge, especially when the approximation of conformal flat- 
ness is eliminated. 



^*We note that for the preferred choice of Wgtatic ^ ~9 the gauge invariant parameters of 
the ISCO found in Ref. [47] do not match well with the same parameter determined via the 
effective-potential method in Ref. [39]. However, if < ^static < lOi the agreement is much 
better. 
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4 Neutron-Star Initial Data 



The construction of initial data for neutron stars requires that the state of the 
neutron-star matter be specified before the gravitational data can be deter- 
mined. Of course, a sohition of the gravitational constraint equations can be 
found in principal for any given energy density p and momentum density j*. 
But, with neutron-star solutions, we are usually interested in situations where 
the matter is in (or nearly in) hydrostatic equilibrium and the gravitational 
fields are also in (or nearly in) equilibrium. 



4.1 Hydrostatic Equilibrium 

For a neutron star to be in true equilibrium, the spacetime must be stationary 
as discussed in §2.4. This means that the spacetime possesses both "temporal" 
and "angular" Killing vectors (cf. Ref. [35]). If the matter is also to be in 
equilibrium, then the 4- velocity of the matter u'^ must be a linear combination 
of these two Killing vectors. If we use coordinates as defined in (53) with the 
angular Killing vector in the (j) direction, then 

m'' = u* [1,0,0, f2]. (113) 

Here, u* and O are functions of r and 9 only. Q. is the angular velocity of the 
matter as measured at infinity. 

It is common to define v as the relative velocity between the matter and a 
normal observer (often called a zero angular momentum observer) so that 

^ = -n^ui^ = au\ (114) 



Vl — V 

The velocity v is then fixed by the normalization condition u^^u^ = — 1. 

If we assume that the matter source is a perfect fluid, then the stress-energy 
tensor is given by 

T^'" = {e + P)uV + Pg^''', (115) 

where e and P are the total energy density and pressure, respectively, as mea- 
sured in the rest frame of the fluid. The vanishing of the divergence of the 
stress-energy tensor yields the equation of hydrostatic equilibrium (often re- 
ferred to as the relativistic Bernoulli equation). In differential form, this is 

dP - (e-hP)(dlnu* - w*W0dr2) = 0. (116) 

If the fluid is barytropic^^, then we can deflne the relativistic enthalpy as 

dP 



h{P) = exp 



e + P 



(117) 



i^For a barytropic fluid, the entropy per baryon and the fractional abundances of the dif- 
ferent nuclear species are determined uniquely by the distribution of baryons. In this case, 
the total energy density e can be expressed as a function of the pressure P (cf. Ref. [35]). 
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and rewrite the relativistic Bernoulli equation as 

ln/i(P) = ln/io + lnM*-ln4- / u^u^dQ. (118) 

J no 

The constants /iq, Mq, and CIq arc the vahics their respective quantities have at 
some reference point, often taken to be the surface of the neutron star at the axis 
of rotation. When uniform rotation is assumed {dfl = 0), Eq. (118) is rather 
easy to solve. The case of differential rotation is somewhat more complicated. 
An integrability condition of (116) requires that u^u^ be expressible as a function 
of fl, so 

= F{n). (119) 

F{^1) is a specifiable function of CI which determines the rotation law that the 
neutron star must obey [35]. 



4.2 Isolated Neutron Steirs 

The simplest models of isolated neutron stars are static (i.e., nonrotating), 
spherically symmetric models that can be constructed, given a suitable equation 
of state, by solving the Oppenheimer-Volkoff (OV) equations [84]: 

_(e + P)(m + 47rr3p) 
r(r — 2m) 

47rr^£ (120) 

a{m + Anr^e) 
r{r — 2m) 

for < r < i?, where R is the radius of the surface of the star. Here, r is an 

areal radius and TO(r) is the mass inside radius r. Exterior to the surface of 
the star, the metric is the standard Schwarzschild metric as in Eq. (59) with 
M = m{R). Interior to the surface of the star, the metric is 

ds^ = -a'df + (l - ^) dr^ + r\de^ + sin^ M</>2). (121) 

The boundary conditions are that m(0) = 0, £(0) is some chosen constant eq, 
and a{R) = y^l — The solutions of this equation form a one-parameter 
family, parameterized by £o which determines how relativistic the system is. 
A method for solving these equations in both the areal coordinate r and an 
isotropic radial coordinate r can be found in Ref. [32] . 

More generally, isolated neutron stars will be rotating. If the neutron stars 
are uniformly rotating, then, for any given equation of state, the solutions form 
a two-parameter family. These models can be parameterized by their central 
density, which determines how relativistic they are, and by the amount of rota- 
tion. If the models are allowed to have differential rotation, then some rotation 
law must be chosen. 



dP 
dr 
dm 
dr 
da 
dr 
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To construct a neutron-star model, the equations for a stationary solution 
of Einstein's equations outlined in §2.4 must be solved self-consistently with the 
equations for hydrostatic equilibrium of the matter outlined above in §4.1. The 
equations that must be solved depend on the form of the metric chosen, and 
numerous formalisms and numerical schemes have been used. An incomplete list 
of references to work on constructing neutron-star models include [102, 22, 34, 
32, 35, 33, 53, 62, 63, 43, 45, 44, 97, 21, 57, 16, 19]. Further review information 
on neutron-star models can be found in Refs. [96, 52, 51]. 

4.3 Neutron-Stcir Binaries 

Neutron-star binaries, like any relativistic binary system, cannot exist in a true 
equilibrium configuration since they must emit gravitational radiation. But, 
as is true for black-hole binary data, for orbits outside the innermost stable 
circular orbit, the gravitational radiation reaction time scale is much longer 
than the orbital period and it is a reasonable approximation to consider the 
stars to be in a quasiequilibrium state. 

A binary configuration obviously lacks the azimuthal symmetry that was 
assumed in the discussions of stationary solutions of Einstein's equations in §2.4 
and hydrostatic equilibrium in §4.1. Fortunately, the condition of hydrostatic 
equilibrium requires only the presence of a single, timelike Killing vector. With 
the assumption that gravitational radiation is negligible, we can assume that 
the matter is in some equilibrium state as viewed from the reference frame that 
is rotating along with the binary. That is, if the binary has a constant orbital 
angular velocity of fi, then the time vector in the rotating frame t''^ is a Killing 
vector^^ and it is related to the time vector in the rest frame of the binary 

by 

t'l' = t'' + n^'', (122) 

where is a generator of rotations about the rotation axis^''. 

Two equilibrium states for the matter have been explored in the literature. 
The simplest case is that of co-rotation, where the 4-velocity of the matter is 
proportional to t'^. In this case, the matter is at rest in the frame of reference 
rotating with the binary system, the corotating reference frame. The second 
equilibrium state is that of counter-rotation, where there is no rotation in the 
rest frame of the binary. We will explore these two cases further below. 

Stationarity of the gravitational field, unlike hydrostatic equilibrium, re- 
quires the presence of separate timelike and azimuthal Killing vectors. For the 
case of a binary system, there is no unique definition of quasiequilibrium. The 
earliest work on constructing quasiequilibrium solutions of Einstein's equations 
stems from work by Wilson and Mathews [103, 104, 105], and others have ex- 
plored similar schemes [11, 18, 12]. Although written in slightly different forms, 
the system of equations for the gravitational fields in all of these schemes are 

^^Bonazzola et al. [18] refer to this as a helicoidal Killing vector. 

^^In spherical coordinates, f = d/d<f>. In Cartesian coordinates, rotation about the z axis 
would be represented by = (0, —y, x, 0). 
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fundamentally identical. While they were developed before the thin-sandwich 
decomposition of §2.3, the thin-sandwich decomposition (see Eq. (51)) offers 
the easiest way to interpret this approach. We consider ourselves to be in the 
corotating reference frame so that our time vector is t''^ . To make the transition 
back to the rest frame of the binary as easy as possible, we write the shift vector 
of our 3+1 decomposition as 

B' = p' + nC, (123) 

so that /?* remains as the shift vector of the 3-1-1 decomposition made with 
respect to the rest frame of the binary system. 

The primary assumptions are that the conformal 3-metric "fij is flat, the 
initial-data slice is maximal so that K = 0, and = 0. We see from (42) that 
the last assumption implies that the conformal 3-geometry is instantaneously 
stationary as seen in the corotating reference frame. The final choice that must 
be made is for the conformally rescaled lapse a. An elliptic equation for the 
lapse can be obtained by demanding the trace of the extrinsic curvature K also 
be instantaneously stationary in the corotating reference frame. This is the so- 
called maximal slicing condition on the lapse. For the particular assumptions 
we have made here, this equation can be written 



(124) 



It is interesting to note that, for a conformally flat 3-geometry, 

{Lby^ = {Lpy^, (125) 

so ft does not appear in the equations for the gravitational fields except in the 
matter terms and possibly in boundary conditions. Equations (51) and (124) 
can be solved for the gravitational fields on an initial-data hypersurface, given 
values for the matter terms and appropriate boundary conditions. 

// wc choose the matter so that it is in hydrostatic equilibrium with respect 
to the pseudo-Killing vector f'^ , then these equations for the gravitational fields 
will yield data that are in quasiequilibrium in the sense that 7^ and K are both 
instantaneously stationary. 

For corotating binaries, the matter is at rest in the corotating reference frame 
of the binary. It is rigidly rotating and hydrostatic equilibrium is specified by 
solving the relativistic Bernoulli equation (118), with dO = 0, self-consistently 
with the equations for the gravitational fields. 

For counterrotating binaries, the matter is not rotating in the rest frame 
of the binary. Counterrotating equilibrium configurations can be obtained by 
assuming the matter to have irrotational flow [98, 91]. As long as the flow is 
isentropic, we can express the enthalpy (117) as^* 

e + P 

h=-^, (126) 
Po 



For isentropic flow, the thermodynamic identity reduces to dh = ^dP. 



35 



where po is the rest-mass density. For irrotational flow, the vorticity of the fluid 
(cf. Ref. [98]) is zero. Combining this with the Euler equation, we find that the 
4-velocity of the fluid can be expressed as 

hui, = Vf,ip, (127) 

where is the velocity potential, or flow field. Equation (127), together with the 
normalization condition u^u^ = — 1, automatically satisfies the Euler equation 
and we are left with the continuity equation which must be satisfied, 

^^^{poun = 0. (128) 

The continuity equation (128), the normalization condition, and Eq. (127) yield 

VV = -(V^¥')V^ Hpa/h) with h = V-(V'^^)(V^(^). (129) 

Stationarity (or quasistationarity) in the corotating reference frame requires 
that 

/lu^i''* = -C, (130) 

where C is a positive constant. Now, in terms of the 3 + 1 decomposition with 
as the shift vector (see Eq. (123)), we find that the Bernoulli equation is 
written 

= -(VV)Vi<^ + B^Vjipf . (131) 

The fiow field <^ must satisfy 

2,,_Riv.YA^_A^ = .^VV- Ar^^ V.ln(^^ 



V^-^'V.^^J--/. = -^VV-^i^^JV.ln^^J (132) 

A = C + B'V,ip 
subject to the boundary condition^^ at the surface of the flow 



0. (133) 



Solving Eqs. (131) and (132) sclf-consistcntly with the equations for the gravi- 
tational fields yields a counterrotating binary in hydrostatic equilibrium. 

As mentioned above, the earliest work on neutron-star binaries was carried 
out by Wilson and Mathews [103, 104]. Wilson et al. [105] describe their ap- 
proach for generating initial data for equilibrium neutron-star binaries. In these 
early works, the equation of hydrostatic equilibrium was not used. Rather, an 
initial guess for the density profile was chosen and the full hydrodynamic system 
was evolved with viscous damping until equilibrium was reached. During each 
step of the hydrodynamic evolution, the equations for the gravitational fields 
were resolved. The resulting data represented neither strictly co- nor counter- 
rotating binary neutron stars. This work led to the controversial result [75] that 

^^This boundary condition comes from the fact that the fluid motion at the surface of the 
star must be tangent to the surface, u'^V^po = 0. 
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each neutron star in the binary may become radiaUy unstable and collapse prior 
to the merger of the pair of stars. While an error was found in this work [50, 76] 
with the result that the signature of collapse is significantly weaker, the contro- 
versy has not yet been completely resolved. 

The first use of corotating hydrostatic equilibrium with the Wilson-Mathews 
approach for specifying the gravitational fields was by Cook et al. [46] for the 
test case of an isolated neutron star. This approach was then used to study 
corotating neutron-star binaries by Baumgarte et al. [11, 12] and by Marronetti 
et al. [71]. Interestingly, turning-point methods for detecting secular instabili- 
ties [94, 95, 54] can be applied to the case of corotating binaries [13]. 

Corotating binary configurations are relatively easy to construct. However, 
it is believed that the viscosity of neutron-star matter is not large enough to 
allow for synchronization of the spin with the orbit [61, 14]. But, if the ini- 
tial spins of the neutron stars are not too large, close binaries should be well 
approximated by irrotational models. Bonazzola et al. [18] (as corrected by 
Asada [7]) developed the first approach for constructing counterrotating binary 
configurations. However, simpler formulations of irrotational flow were devel- 
oped independently by Teukolsky [98] and Shibata [91], and Gourgoulhon [56] 
showed that all three approaches were equivalent. Numerical solutions of the 
equations for irrotational flow coupled to the equations for the gravitational 
fields are more difficult to construct than those for corotation because of the 
boundary condition (133) on the flow field that must be applied on the surface of 
each neutron star. This boundary condition is particularly difficult to implement 
because the location of the surface of the star is not known a priori, and will 
move as the equations are being solved. The first models of irrotational binary 
neutron stars were constructed by Bonazzola et al. [20], Marronetti et al. [72], 
and Uryu and Eriguchi [101]. A full description of the numerical methods used 
by Bonazzola can be found in Ref. [58]. 
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